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Abstract 

Using data assimilation framework, to merge information from model 
and measurement, an optimal reconstruction of the neutronic activity 
field can be determined for a nuclear reactor core. In this paper, we 
focus on solving the inverse problem of determining an optimal repartition 
of the measuring instruments within the core, to get the best possible 
results from the data assimilation reconstruction procedure. The position 
optimisation is realised using Simulated Annealing algorithm, based on the 
Metropolis-Hastings one. Moreover, in order to address the optimisation 
computing challenge, algebraic improvements of data assimilation have 
been developed and are presented here. 
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1 Introduction 

Data assimilation methodology allows to build optimal reconstruction of activ- 
ity field, within a nuclear core, using information coming from both model and 
measurements [TJ HI |3J- The efficiency of data assimilation for physical field 
reconstruction has already been demonstrated in several articles in meteorology 
[HOE]. Demonstration of its efficiency for nuclear core activity field recon- 
struction was also done in recent articles such as pj |H1 H] ■ The main points of 
data assimilation are presented in appendix |B) 

However, in all those applications of data assimilation, the measurement 
network is considered to be known. The inverse problem consists then in op- 
timizing the location of the measuring instruments, in order to get the best 
possible results from the data assimilation reconstruction procedure. Unlike the 
meteorological domain, for which the question cannot be addressed due to the 
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great size of the problem (over 10 6 measurement locations), the issue can be ad- 
dressed and resolved for the nuclear core, as only about 10 2 to 10 3 measurement 
locations have to be taken into account. 

The chosen method to optimise the instrument location is the well known 
Simulated Annealing algorithm, based on a Metropolis-Hastings one [TUJ [TT] . 
This algorithm was originally designed to deal with problems of statistical 
physics. The first version, from N. Metropolis [TP] , was only focused on Boltz- 
mann equation. Then W. K. Hastings generalized it to other cases [TT]. From 
these methods, S. Kirkpatrick proposed some improvements [12] by changing the 
definition of the energy, used in the algorithm to optimise position in system 
presenting an inner local order. Other applications of this algorithm already 
exist for example in nuclear incident network design since a few years |131 1 14] . 
A brief description of the main points of Simulated Annealing is presented in 
appendix [A] 

One of the key point in the Metropolis-Hastings algorithm is the permutation 
of an instrumented location in the core with a non instrumented one. This can 
be numerically costly to redo data assimilation calculations for each case, as it 
includes big matrix inversions. To make those calculation feasible, we developed 
methods of fast calculation, that treat instrument movement as a sequence of 
loss and gain of instruments. This can be considered as a perturbation of the 
initial state, which become then rather cheap to compute. The instrument loss 
fast calculation was already used in [5]. The algebraic acceleration method is 
recalled in appendix [D] 

In the paper, the first section describes the general framework of the study. 
In a second section, the parametrisation of simulated annealing method is pro- 
vided. Then the optimisation results, of simulated annealing coupled with an 
assimilation procedure, are shown, either starting from standard reactor in- 
struments distribution or from a random distribution. In a next section, the 
synthesis of both result is done. The global conclusion on the result of instru- 
ment location optimisation using simultaneously simulated annealing and data 
assimilation is then given. 

2 Definition of the nuclear core working frame- 
work 

The experimental set up of this study is the standard configuration of a 900 
MWe nuclear Pressurized Water Reactor (PWR900), and we study the neu- 
tronic activity field reconstruction. The core is filled with a total of 157 vertical 
nuclear fuel assemblies, among those 50 are instrumented with Mobiles Fissions 
Chambers (MFC) to measure the neutronic activity fields. An horizontal slice 
of a PWR900 core is represented on the Figure [TJ The question is then to find 
the location of the 50 instruments within the 157 allowed positions. We don't 
consider here any constraint on the location of measurements, despite the fact 
that real constraints exists (mainly, MFC and control rods, used for safety and 
steering, are incompatible). It is not realistic, but, without these constraints, 
the optimisation problem is more similar to experimental capabilities of new in- 
struments on latest reactors, and is far larger and difficult to solve. Solving this 
unconstrained optimisation problem is then a prerequisite to take into account 
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any constraint. 
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Figure 1: In the nuclear core, MFC instruments position is indicated by assem- 
blies in black within an horizontal slice of the core. The assemblies without 
instrument are marked in white, and the neutronic reflector, out of the reactive 
core, is in grey. 

To perform data assimilation, both simulation code and experimentally mea- 
sured data are needed. For the simulation, the EDF newest experimental cal- 
culation code COCAGNE for nuclear cores [T5J [TB] is used in a standard con- 
figuration. A description of basic features of the underlying neutronic model is 
done in the appendix [Cj The physical assemblies are numerically represented 
using 29 vertical levels. Thus, the size of the data assimilation control vector x 
is 4553 (157 x 29). The size of the observation vector y° is 1450 (50 x 29). All 
the details of data assimilation parameters, derived from those choice of model 
and observation space, are given in appendix [C] 

To have a good understanding of the instrumentation effect, we study various 
scenario of instrument configurations (even some that do not exist operationally 
and so cannot be tested experimentally). For that, synthetic data are used, that 
allows to have an homogeneous approach all along the document and to use twin 
experiment methodology. Synthetic data is generated from a model simulation, 
filtered through an instrument model, and noised according to a predefined 
measurement error density function (usually of Gaussian type) . 

3 Parametrisation of the simulated annealing al- 
gorithm 

The global description of the well known simulated annealing and the Metropolis- 
Hastings algorithms is done in the appendix [X] This section only focus on the 
specific parameters related to the problem of instrument locations optimisation, 
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coupled with a data assimilation procedure for the activity field reconstruction. 

3.1 Definition of the energy 

In the present case, the aim is not to minimise a real energy, but to improve the 
quality of reconstruction x a of the neutronic activity fields on the nuclear core. 
The reconstructed field x a is the result of a data assimilation procedure, using 
parameters described in appendix [C] for a given instrumental configuration. 

Thus we would like to minimise the distance to the true value x* of the 
core activity fields, that we know because we decide to set ourselves in twin 
experiment framework. The energy E will then be defined as the norm of the 
difference between the analysis x a and the true value x*, written as the following 
formula: 

£(x a ) = ||x a -x*|| (1) 

This definition of the energy will be used, all along the Metropolis-Hastings 
algorithm, for the required term in formula [7| 

3.2 Measuring the evolution of the quality 

To measure the evolution of the quality of an instrument configuration, as a 
function of the iteration in the process, a natural choice is to plot the value of 
the energy E(x. a ). 

However, this value as no real physical meaning. Thus, it is more interesting 
to study the difference of i?(x a ) with respect to a chosen reference E(x£ e f). This 
is only a shift in the energy function, but this leads to a better suited function 
to represent the quality of the x a reconstruction: the best is the instruments 
location, the closest the new quality is to zero. 

We choose this reference value to be given by an analysis x" e *, obtained 
by data assimilation, assuming each assembly in the core is instrumented with 
one MFC. In a idealized case, such a "limit state" can be accessed by experi- 
ment. The calculation is done with respect to the same parametrisation of data 
assimilation as described in the appendix [C] 

The quality evaluation factor q is then given by the following formula: 

g(x a ) = £(x a ) - E^ ref ) = ||x Q - x'|| - |K e/ - x*|| (2) 

Its values are negative. Thus, if the Metropolis-Hastings algorithm improves 
the instrument location choices, we expect to get an increasing curve for q over 
the iterations. 

3.3 Starting state for the optimisation algorithm 

The Metropolis-Hastings algorithm is an iterative algorithm, requiring a starting 
point or instrument state to begin optimisation. We choose particular starting 
points to illustrate how is working the algorithm. Moreover, the Metropolis- 
Hastings algorithm is also a stochastic one. Thus assuming no seed is chosen 
in the random number generator, two realisations of the algorithm will lead to 
different sequences of evolution. 

The classical PWR900 configuration is expected to be a very good one, re- 
liable and being designed years ago. Thus, as a first possibility for the starting 
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point of the Metropolis-Hastings algorithm, this configuration can be chosen. 
This will assess the quality of this design (that originally satisfy also to a lot of 
operational constraints, on the contrary to the one we are seeking by optimiza- 
tion). In order to illustrate how is working the algorithm, two different random 
realisations of the algorithm will be presented based on the same initialisation 
in section |4j 

A second possibility for the starting point is to randomly choose it. This will 
allow comparison of the optimal search with the one obtained with PWR900 
configuration as initial point. 

An another comparison can be done by using less instruments than the 
50 classically considered, in order to look at the influence of the number of 
instruments. In this case, the starting point has to be randomly chosen, as 
there is no standard repartition of instruments. 



3.4 Other parameters of the algorithm 

As is stated in appendix [Aj parameters need to be chosen to run a Metropolis- 
Hastings algorithm. There are mainly 2 parameters to set in the present case. 

The first one is the maximal number of iterations, which is the number of 
time that we go through the loop on instrumentation swapping. This number 
will be set to 1800. This value was chosen for two reasons. The first point is that 
we notice, through the setting up trials, that with such a number of iterations, 
we reach a steady state in every case we processed. As a second point, it 
corresponds to reasonable computing time of about 10 hours (depending on the 
computer used) per search and per processor. 

The second parameter is the choice of the thermalisation function that fixes 
the evolution of the pseudo-temperature in the algorithm. We experiment sev- 
eral standard thermalisation formulation, and put our final choice on a classical 
simple inverse function of the iteration, defined as follows: 

Ti = ^r (3) 
i + 1 

where i is the iteration step. The initial value To of the temperature was set 
to 0.05. With such a definition, we got satisfactory result for the optimisation 
algorithm. 



4 Result of optimisation from PWR900 config- 
uration 

In the figure [2] and [3] two cases are presented, based on the same initialisation 
and on different random realization. On those figure, the red crosses represent 
the quality of the configuration at the current iteration, and the blue asterisks 
represent the quality of the best state founded since the beginning. 

First, we are looking at the best state found curve. In both figures [2] and 
[3j it can be noticed that the algorithm reach a better state at the end than 
the initial one. Looking at the quality of the final state in the two cases, it 
can be noticed that they are very close. Looking on a set of realisation of the 
optimisation (not presented here), the final state is always around the same 
level of quality with a rather small dispersion (around 0.002, that is around 
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Figure 2: Evolution of the quality of the reconstruction as a function of the 
iteration on 1800 steps (first random realization). At each step are represented 
the best evaluated quality as well as the quality at the current iteration. 

Optimisation from PWR900 
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Figure 3: Evolution of the quality of the reconstruction as a function of the it- 
eration on 1800 steps (second random realization). At each step are represented 
the best evaluated quality as well as the quality at the current iteration. 



5%). In both cases, in a first phase, the quality of the configuration improves 
fairly quickly and a good configuration can be found after around 300 iterations. 
After that, the improvement of the configuration is very slow and even stagnate 
above 800 iterations. 

Then, it is interesting to have a look at the value of the quality of the current 
state for each iteration, to better understand the behaviour of the Metropolis- 
Hastings algorithm. The plot of the current state in figure [2] is the typical one we 
got over several realisation. The current state is " oscillating" around the current 
best state, to find a new optimal state. The behaviour of the algorithm plotted 
in figure[3]is less common, but very illustrative of the potentiality of the method. 
We notice that the quality of the final state is about the same as the one of the 
figure [2j But in figure [3j we notice also that the current state can get very far 
away of the best state, for a very long time steps. This kind of minimum search, 
far from the current state, is one of the strong point of the present algorithm. 
This research of a new optimal state is possible because of the highest pseudo 
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temperature at the beginning of the process given by equation [3j The most 
interesting fact is that, whatever how far the method explore the state space, 
finally it converge again to a state of good quality. This property of Metropolis- 
Hastings algorithm, noticed qualitatively here, have been demonstrated rather 
recently mathematically as mentioned in reference [17] . 

At the end of the optimisation process, a new set of locations for the instru- 
ments is obtained. The set of location coming from the two previous optimisa- 
tions are plotted in figures [3] and [5] 

Location set represented inland [5] have, at glance, no common point. How- 
ever, from the point of view of the quality defined by equation [2] those two 
repartitions are very close. This prove not only that one optimal state does ex- 
ist, but also that there is an ensemble of optimal states of rather close quality. 

Examining more precisely both figure [4] and [5] (and other realisation not 
presented here), some common features can be discussed. In particular, there is 
a tendency to form small clusters of instruments, mainly localised in the centre 
of core. Another feature is the repartition of several individual instruments 
around the core. However, a detailed study over many realisations shows that 
such features are not significant. 

5 Result of optimisation from random configu- 
ration and dependency to instruments num- 
ber 

Without an initial guess of the instrument repartition as the one of standard 
PWR900, the natural choice of a starting repartition is to take them randomly. 
To be in the same condition as in the PWR900 case, we choose to use 50 
instruments randomly located in the core. Using random location is, at the 
same time, a good opportunity to get information on method and result when 
less instruments are available. Thus the case where only 10 instruments are in 
the core will be studied too. 

To define the problem, a brief study of the property of the result of data 
assimilation quality for random instruments repartition has been done. Over a 
set of 1000 random distribution of the 50 instruments, it can be noticed that the 
quality q of all data assimilation reconstructions x a (with respect to formula [2| 
are distributed following as Gaussian distribution of mean m — —0.0832 and a 
standard deviation a — 0.0078. 

This already means that, with a value of q = —0.0476, the PWR900 standard 
repartition is an excellent choice, with a far greater quality. 

Assuming only 10 instruments, the distribution of the quality q of the anal- 
ysis x a still have Gaussian shapes but with a mean value of m — 0.0868 and a 
standard deviation of a = 0.076. The very small difference between the mean 
value of both distributions proves that, due to distance of influence included in 
data assimilation procedure, inappropriate localisation of the instruments can 
result in decreasing performance of the method, as already noticed in reference 

Now we look at the evolution of the quality of the data assimilation recon- 
struction as a function of iterations of Metropolis algorithm. In the figures [6] and 
[7j those two cases are presented. Within those figures, as in figures [2] and |3j the 
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Figure 4: The optimised positions of MFC instruments after a first random 
realisation are localised in assemblies in black, within an horizontal slice of the 
core. The assemblies without instrument are marked in white and the reflector, 
out of the reactive core, is marked in grey. 




x position 

Figure 5: The optimised positions of MFC instruments after a second random 
realisation are localised in assemblies in black, within an horizontal slice of the 
core. The assemblies without instrument are marked in white and the reflector, 
out of the reactive core, is marked in grey. 



crosses in red represent the quality of the configuration at the current iteration 
and the asterisks in blue represent the quality of the best state founded. 

In both figures [6] and [7J it can be noticed a very fast improvement, of the 
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Figure 6: Evolution of the quality of the reconstruction using 50 instruments, 
initially randomly localized, as a function of the iteration on 1800 steps. At 
each step are represented the best evaluated quality, as well as the quality at 
the current iteration. 

Random Start 10 instruments 

-0.01 

-0.02 

-0.03 

-0.04 

| -0.05 
O 

-0.06 
-0.07 
-0.08 
-0.09 

200 400 600 800 1 000 1 200 1 400 1 600 1 800 

Iteration 

Figure 7: Evolution of the quality of the reconstruction using 10 instruments, 
initially randomly localized, as a function of the iteration on 1800 steps. At 
each step are represented the best evaluated quality, as well as the quality at 
the current iteration. 



quality of the instruments setting, in very few iterations of the Metropolis- 
Hastings algorithm. This initial phase is followed by a slow evolution of the 
quality, and then by a stagnation phase. Stagnation phase is reached around 
after 600 iterations (resp. 300) for 50 (resp. 10) instruments. This difference 
can be easily explained, considering the space of possible configurations for 
50 instruments in the core is far more large than the one of 10 instruments. 
The final state quality value is of -0.0305 (resp. -0.0504) for 50 (resp. 10) 
instruments. Those qualities after optimisation are coherent with the amount 
of available information, that is, of number of available instruments. This was 
shown in [S]. 

The behaviour of the current state is typical of the Metropolis-Hastings 
algorithm on both figures [6] and [7| It is then interesting to compare the locations 
of the instruments in the optimised state for both cases of 50 and 10 instruments. 





9 



Those locations are presented in figures [8] and [9] 
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Figure 8: The optimised locations of 50 MFC instruments are localised in as- 
semblies in black within the horizontal slice of the core. The assemblies without 
instrument are marked in white and the reflector, out of the reactive core, is in 
grey. 




2 4 6 8 10 12 14 

x position 



Figure 9: The optimised locations of 10 MFC instruments are localised in as- 
semblies in black within the horizontal slice of the core. The assemblies without 
instrument are marked in white and the reflector, out of the reactive core, is in 
grey. 

In both figures [8] and [9] appear the same features as the one obtained in 
figures [4] and [5j coming from optimisation of an initial PWR900 standard con- 
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figuration. Those features are a tendency to form small cluster of instruments 
mainly localised in the centre of core, and the repartition of several instrument 
around the core. However, this is still statistically irrelevant. 



6 Synthesis of the results 



Using the information gathered from the various uses of the method, we make 



a synthetic plot of the information obtained and comment it. On figure 10 



are 



plotted the quality of the best evaluation as a function of the number of iteration 
coming from figures [2] [6] and [7j 



Synthesis of results 



50 instruments PWR900 
50 instruments random 
10 instruments random 




Figure 10: Evolution of the quality of reconstruction as a function of the it- 
eration on 1800 steps, with 50 instruments initialised by PWR900 or random 
configurations, and with 10 instruments initialised randomly. 



The first remark on figure 10 is that, whatever are the initial condition 
(instruments repartitions or number of instruments), the simulated annealing 
method allows to find an better solution. 

One interesting result of this study is, as expected, the good quality of the 
PWR900 standard instruments distribution. Without any optimisation, the 
quality of this repartition is q = —0.0476, that should be compared to a mean 
value of m = —0.0832 for several random repartition with 50 instruments. We 
recall that the optimisation procedure is done without any constraints, which 
is not physically comparable with the PWR900 and its standard instruments 
distribution. 

Thus, this is not a good idea to chose randomly the position of the instru- 
ments, even with a data assimilation procedure to make an optimal determina- 
tion of the activity field. However, looking for the case with 50 instruments, 
after optimisation, the final results have the same quality with a quality factor 
q close to —0.0325. This is a experimental confirmation of the theoretical result 
about convergence of simulated annealing method |17j . this even if the final 
state for the same quality do not look alike as it can be seen in figures [3] and [8] 

Using different number of instruments as in section |5j we notice interesting 
results. Observing the starting quality value q of results with data assimilation, 
the situation is rather blurry, as the quality with 50 instruments is close to 
the one with 10 instruments. Then, after position optimisation by simulated 
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annealing, ordering is clearly done according to the instruments number: the 
more instruments we use, the best quality we get. Moreover, this result is 
general, and final quality is always ordered respect to the number of available 
instruments. 

All those results demonstrate the efficiency and robustness of the simulated 
annealing method for instruments position optimisation. This is even more 
stringent considering that the method is applied in link with data assimilation 
technique, that already enrich the global information on the whole system. 

7 Conclusion 

From this study, two important results are enlightened. 

First, the standard PWR900 instruments repartition is characterized by an 
excitingly good quality. This instruments repartition is a priori the best we 
know, even if it was not originally designed using a data assimilation framework 
background. 

The second point is that the simulated annealing method can always find a 
ameliorated instrument locations set. Moreover, as expected theoretically [T7] . 
the final state is of the same quality whatever is the starting point. And using a 
reduced number of instruments, the final organisation has always a quality that 
is increasing with the increase of the number of instruments. 

Those results demonstrate that, within the framework of neutronic simula- 
tion for nuclear PWR, and using a reasonable size for the discrete space in which 
to optimize the instruments (below 10 4 points), it is possible to do instruments 
repartition optimisation with data assimilation method. Such a work can be 
done more efficiently considering an algebraic optimisation of the calculations, 
and become still very reasonable in computing time. 

Within the framework of nuclear core neutronic simulation, several devel- 
opments can be done from those results. Specially, instruments repartition can 
be done under constraints such as symmetry or specific position for example. 
Moreover, use of data assimilation open the way to obtain simulated annealing 
optimisation with heterogeneous available instruments. 
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A Metropolis-Hastings algorithm 



The algorithm used for position optimisation is the Metropolis-Hastings one. 
It was initially used to describe the thermodynamically evolution of a system, 
and the first algorithm was written by N. Metropolis, A. W. Rosenbluth, M. N. 
Rosenbluth, A. H. Teller, and E. Teller [10] • In this version, only a Boltzmann 
type distribution was considered as it is the main one used in statistical physic. 
Then, in 1970, W. K. Hastings extends this algorithm to other distributions 

P. 

The method is a random walk in instruments localisation space, and so an 
iterative algorithm. This random walk use a probability Q to chose the next 
possible position x, moving from x^ at step k. The transition probability is 
Q(xk,x). Thus the probability that the whole system go from Xk to x at step 
k + 1 (which means Xk+i — x) is then defined as follow: 

P(x t+1 = x\x t ) =mm\ v -,U (4) 

[■K(Xt)Q{X,Xt) J 

The remaining probability is the one to stay in Xk- 

P{x t +\ = x t \x t ) = 1 - P(x t+ i = x\x t ). (5) 

In the very common case of symmetric evolution, which means that Q(x, y) = 
Q(y,x), some extra simplifications can be obtained. When optimising instru- 
ments location, we are in such changing position of an instrument 
and then making the invert operation leads to the same first repartition. More- 
over, the choice of both, the instrument to move, and its new localisation, are 
perfectly random: 

• with a probability of 1 if tt(x) > ir(xt) 

• with a probability of else 

using a given probability distribution tt. Any kind of probability distribution can 
be taken into account. In the present case, like in many applications of Metropo- 
lis algorithm, a Boltzmann distribution will be used as probability function. The 
Boltzmann distribution is defined by the following formula: 

-E(x) 

n(x) = e t (6) 

where E(x) = E is an evaluation of the quality of the state (such as a physical 
energy or something else) of the state x, and where T is a temperature or a 
pseudo-temperature associated to E. Thus, in the equation [4j when we make 
the ratio , we obtain the following formula: 

7r(x t ) ' o 

-E(x)+E(x,.) -A(x,x fc 

7r(x) = e t = e t (7) 

This is this formula used here for the Metropolis algorithm. 

The final algorithm of positions optimisation is then a loop, that last until 
a condition on number of iteration or on energy limit is reached. Within the 
loop, the following actions are realised: 
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• Swap two instruments or displace one instrument by swapping with a void 
place, 

• Evaluate the pseudo energy E of the new configuration, 

• Compare and record the best configuration founded, 

• Chose evolution of the instrument pattern according to probability law 
given by equation [7J 

• Go back to first step. 

The algorithm is finally simple and reliable. Those type of methods have 
been used for positions optimisation in several applications, and one of the first 
and most relevant result was obtained for optimisation of chip design by S. 
Kirkpatrick, C. D. Gelatt and M. P. Vecchi [12]. 

B Data assimilation 

We briefly introduce the useful data assimilation key points to understand their 
use, as applied in p i [T9 l 120 ] . 

Data assimilation is a wide domain and these techniques are, for exam- 
ple, the keys of the nowadays meteorological operational forecasts [21] . This is 
through advanced data assimilation methods that weather forecasting has been 
drastically improved during the last 30 years. All the available data, from satel- 
lites, aircraft or ground measurements, are used in conjunction with complex 
numerical weather models. 

The goal of data assimilation methods is to estimate the true value x* of the 
state of the considered system, where the t index stands for " true" . The basic 
idea is to combine information from an a priori knowledge on the state of the 
system (usually denoted as x b , with b for "background"), and some measure- 
ments (referenced as y°, with o for "observations"). The background is usually 
the result of numerical simulations, but can also be derived from any a priori 
knowledge. The result of data assimilation is called the analysis, denoted by 
x a , and it is an estimation of the true state x* we look for. 

The control and observation spaces are not necessarily the same, and a 
bridge between them has to be built. This is the observation operator H, that 
transforms values from the space of the background to the space of observations. 
For data assimilation purpose, we use the linearised operator H of H around 
the background x^. The reverse operator, converting observation increments to 
background increments, is given by the transpose H T operator of H. 

Two other ingredients are necessary. The first one is the covariance matrix R 
of observation errors, defined as R = S[(y°-i/(x*)).(y°-iJ(x*)) T ], where E[.} 
is the mathematical expectation. It is obtained from the known errors assuming 
unbiased measurements, which means E[y° — H(x t )] = 0. The second one is 
the covariance matrix B of background errors, defined as B = E[(x b — x 1 ).^ — 
x*) T ]. These errors on the a priori state are also assumed it to be unbiased. 
There are many ways to get this a priori state and background error matrices. 
However, those matrices are commonly obtained from the output of a model by 
an evaluation of accuracy, or are the result of expert knowledge. 
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Within this formalism, under a static assumption of state equations, the 
analysis x a is the Best Linear Unbiased Estimator (BLUE), and is given by the 
following equation: 

x a =x b + K(y° - Hx b ), (8) 
where K is the gain matrix such as: 

K = BH T (HBH T + R) _1 . (9) 

Moreover, we can express the analysis error covariance matrix A characterising 
the analysis errors x a — x . This matrix is derived from K as: 

A = (I-KH)B, (10) 

where I is the identity matrix. 

It is worth noting that solving Equation [8] is, if the probability distribution 
is Gaussian, equivalent to minimise the following function J(x), x a being the 
optimal solution: 

J(x) = (x - x 6 ) T B- x (x - x") + (y° - HxjV 1 (y° - Hx). (11) 
This minimisation is known in data assimilation as 3D-Var methodology 18J. 

C Data assimilation implementation 

C.l Brief description of the nuclear core modelling 

The aim of a neutronic code, like COCAGNE [HI [16] from EDF, is to evaluate 
the neutronic activity field and all associated values within the nuclear core. 
This field depends on the position in the core and on the neutron energy. To do 
such an evaluation, the population of neutrons are divided in several groups of 
energy. Classically, two energy groups are considered, describing the neutronic 
flux by $ = ($i,$2)- The material properties depend on the position in the 
core. The neutronic flux $ is identified by solving two-group diffusion equations 
described by: 

' -div(^ lg rad$ 1 ) + (S ol + E r )$i 

= -^iS/i^x + i/ 2 £ /2 $ 2 ) (12) 
t -div(L> 2 grad$ 2 ) + £ a2 $ 2 - S r $i = 

where k is the effective neutron multiplication factor, all the quantities and the 
derivatives (except k) depend on the position in the core, 1 and 2 are the energy 
group indexes, E r is the scattering cross section from group 1 to group 2, and 
for each group, E a is the absorption cross section, D is the diffusion coefficient, 
and vYj f is the corrected fission cross section. 

The cross sections also depend implicitly on the concentration of boron, 
which is a substance added in the water used for the primary circuit to con- 
trol the neutronic fission reaction. This control is described through a feedback 
supplementary model. This model takes into account the temperature of the 
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materials and of the neutron moderator, given by external thermal and thermo- 
hydraulic models. A detailed description of the core physic and numerical solv- 
ing can be found in reference |22j . 

The overall numerical resolution consists in searching for boron concentra- 
tion such that the eigenvalue k is equal to 1, which means that the neutron 
production in the core is stable and self-sustaining. It is named " critical boron" 
concentration computation. 

The activity in the core is obtained through a combination of the fluxes 
<f> = (<f>i, $2). Numerically, the activity is given on a chosen mesh of the core. 
Using homogeneous materials for each assembly (for example 157 in a classical 
EDF PWR900 reactor), and choosing a vertical mesh compatible with the core 
(usually 29 vertical levels), this result is an activity field of size 157 x 29 = 4553 
that cover all the core. 



C.2 The observation operator H 

The observation operator H is, in the present application, a selection proce- 
dure. This procedure extracts the values corresponding to effective measure- 
ments among the values of the numerical model space. The normalisation pro- 
cedure is a scaling of the value with respect to the geometry and power of the 
core. The entire process is linear. Then the linear operator H is identical to H . 
Size of the operator depend directly on the number of instruments available. In 
a PWR900 as we treat it, the size is (1450 x 4553). 



C.3 The background error covariance matrix B 

The B matrix represents the covariance between the spatial errors for the back- 
ground. In order to get those, we estimate them as the product of a correlation 
matrix C by a normalisation factor. 

The correlation C matrix is built using a positive function that defines the 
correlations between instruments with respect to a pseudo-distance in model 
space. Positive functions have the property (via Bochner theorem) to build a 
symmetric defined positive matrix when they are used as matrix generator [231 
|2~4"] . In the present case, Second Order Auto-Regressive (SOAR) function is used 
to prescribe the C matrix. In such a function, the amount of correlation depends 
from the Euclidean distance between spatial points. The lengths of radial and 
vertical correlation (denoted as L r and L z respectively, and associated to the 
radial r coordinate and the vertical z coordinate) have different values, which 
means we are dealing with a global pseudo euclidean distance. The function 
can be expressed as follow: 



The matrix obtained by the above Equation 13 is a correlation matrix. It is 
then multiplied by a suitable variance coefficient to get a covariance matrix. 
This coefficient is obtained by statistical study of differences between model 
and measurements on real case. In our case, the size of the B matrix is related 
to the size of model space, so it is (4553 x 4553). 
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C.4 The observation error covariance matrix R 

The observation error covariance matrix R is approximated by a diagonal ma- 
trix. This assumes that no significant correlation exists between the measure- 
ment errors of the measure instruments. On the diagonal, the usual modelling 
is to take each value as a percentage of the corresponding observation. This can 
be expressed as: 

R.ij = W)i) 2 . vj (u) 

The parameter a is fixed according to the accuracy of the measurement and of 
the representative error associated to the instrument. The size of the R matrix 
is related to the size of observation space, so it is here of (1450 x 1450) for a 
PWR900. 



D Fast calculation method for data assimilation 
on perturbed instruments network 

The direct calculation of influence of the displacement of an instrument, which 
is mandatory for the Metropolis algorithm, is very time consuming because of 
the size of the matrices. Thus, to shorten this computing time, we developed 
special methods to calculate the removal or the addition of an instrument as a 
perturbation of a known state. Thus, changing the position of an instruments 
in the network corresponds to remove an instrument and then to add one, or 
inversely. Calculating the new analysis using such approach is then substantially 
faster. We present in this section the two methods, for fast analysis calculation 
when adding or removing instrument. Both techniques are based on Schur 
complement [25] . 



D.l Instrument removal 

Within the BLUE assimilation method, the limiting factor in calculation time 
are the matrices inversions. In Equation |9j the costly part is the inversion of 
the term M defined as: 

M = HBH T + R. (15) 

In the present case, the number of terms of M is around 2.10 6 , so the inversion at 
each step is rather time consuming when iterating for the Metropolis algorithm. 
Then we try to optimise the computing cost. 

We noticed that the calculations are more time consuming when only few 
instruments are removed. In this case the M matrix remains still huge. 

Thus, the idea is to use the information obtained in the inversion of the full 
size matrix to shorten calculation, by calculating smaller size matrix. In this 
case, we want to calculate the new matrix as a perturbation of the original one. 
This is done by exploiting the Schur complement of the matrix. 

We want to suppress some instruments to a given physical configuration. 
With respect to the Equation [9j we need to calculate a new matrix K„. The n 
index is standing to designate the new matrix we want to calculate. According 



to Equation 15 we have to determine a new matrix M„. 

This matrix M n is obtained from the knowledge of the invert of the matrix 
M s calculated over all the instruments. The subscript g is used to designate 
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the reference matrix we start from, according to Equation [l5j in the case of a 
complete initial instruments repartition. 

All the components of the new matrix K„ can be obtained by suppressing 
the lines and columns corresponding to removed instruments in M g , inverting it 
and then multiplying this matrix by the corresponding H„ and B n . We notice 
that, in our case, B„ = B g as procedure do not affect the model space. 

To make the explanation easier (and without loosing generality) , we assume 
that the suppressed instruments correspond to the lower square of M g . If it is 
not the case, it is always possible to reorganise the matrix in such a way. 

Now we put the matrix M g under a convenient form, separating remaining 
measures from removed ones. Assuming the starting matrix M g is of size m x to, 
and assuming we plan to keep p measurements and to suppress s, M g can be 
written in the following way: 

where: 

• P g contains the remaining measurements, and is a p x p matrix, 

• S g contains the suppressed measurement, and is a s x s matrix, 

• Q g ct R g represent the dependence between remaining measured and sup- 
pressed ones. In the particular case we are dealing with, = R g . How- 
ever, no further use of this property is done. 

With such a decomposition, we got the equality to = p + s. The P g matrix 
corresponds to the remaining instruments, thus we have the equality: 

P 3 = M„, (17) 

The Equation [16] gives the decomposition required to build the Schur comple- 
ment of this matrix |25j . Under the condition that P g can be inverted, the 
Schur complement is the following quantity: 

Sg - RgP-'Qg, (18) 

and is noted (M 9 /P s ). This notation reads as "Schur complement of M g by 
P " 

Thus we look for a cheap way to calculate P g 1 knowing M g 1 . For that, we 
use the Banachiewicz formula [5S], that gives invert of M g as a function of P g , 
Q g , R g , S g and (M g /P g ) matrices: 

^ = (rI tV (i9) 

P" 1 + P- 1 Q g (M g /P g )- 1 R g P g 1 -P-iQ ff (M g /P s )- : 

We define the 4 sub-matrices P g , Q g , R g and S g by: 

P g = P g - 1 +P g - 1 Q g (M g /P g )- 1 R g P g - 1 , (20) 

Q g = -P g - 1 Q g (M g /P g )~ 1 , (21) 

R g = -(M g /P g )- 1 R g P g ~ 1 , (22) 

Sg = (M ff /P g ) _1 . (23) 
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Rearranging those terms we get: 

P,- 1 = P 9 - Q g S~ l K g . ( 24 ) 

By hypothesis, the inverse M^ 1 of the global matrix is known, and it is possible 
to extract P s , Q g , R 5 and S g from the whole inverted matrix. Thus, the main 
cost to obtain the inverse of P s , of size p x p, becomes the inversion of the 
matrix S g , which size is q x q. If the number of measurements to suppress is 
smaller than the number of remaining ones, this methods gives a notable gain. 
As soon as the matrix P g 1 = M^ 1 , final calculation of K„ is straightforward. 



D.2 Fast calculation of the gain matrix K in the case of 
instruments gain 

The next step is to calculate the required matrix in the case of instruments gain. 
In this case the situation is more tricky, as the dependency between the various 
instruments need to be taken into account, if they do exist. As previously, the 
technique used is based on Schur complement. 

In this part, we consider a system with two instruments labelled with number 
1 and 2, and characterised by observation operators Hi et H 2 . We note di and 
di the size of the spaces respectively associated with Hi and H2. The size of 
the data space is denoted n. The matrix Hi and H2 are of size d\ x n and 
rf 2 x n respectively. 

The two error covariance matrices on measurement Ri and R 2 are of size 
d\ x di and d2 x di respectively. Without loosing generality, the B matrix can 
be assumed to be unique. Both operators Hi and H2 are application from the 
same space of data toward different observation spaces. The dimension of the 
observation space is n, then the size of this matrix B is n x n. 

Starting from the formula [9] given for independent calculation of the analysis, 
we got : 

Ki =BHf(HiBHf + Ri) 1 

(25) 

K 2 = BH 2 n (H 2 BHj + R 2 ) 1 

The coupled observation operator of the whole measurement system is de- 
fined by combination of the instruments 1 and 2. Thus, the complete observation 
operator can be written as follow: 

H =(5;) <*> 

This is completely equivalent to exchange the two observation operators Hi 
and H 2 . The final matrix obtained is of size (d± + d^) X n. Looking now at the 
structure of the covariance error matrix R, it can be written as follow: 

In the above formula, the R12 sub-matrix is the correlation error matrix between 
the two instruments set. This sub-matrix contains the information required to 
combine both assimilation. Even if we restart a data assimilation from the 
beginning, this information need to be added. Without any information, in 
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many case, it can be assumed that no correlation does exist between the errors 
of each set of instruments. This is equivalent to take R12 = 0. However, to keep 
generality, we will treat the problem assuming a non zero R12 matrix. 

From that point, all the mandatory element to realise data assimilation 
calculation are available. The complete matrix K, described in Equation]!)} can 
be then decomposed according to the above information as follow: 

K = B ( Hf ) x 

*>r ,28> 

Doing the multiplication of Hi and H 2 matrix with B matrix, the previous 
expression becomes: 



K = B ( Hf H 



R^ R 12 \ + / HxBHf HiBH| ^ 1 (29) 
R 12 R 2 / \ H 2 BH 1 H 2 BH 2 



In the term that needs to be inverted, the matrix is decomposed as a sum 
of diagonal element and extra diagonal ones. The diagonals terms correspond 
to previously known elements, and the extra diagonal ones appear due to the 
interdependency between variables, both in the model (with terms HiBH^ and 



H 2 BH^) and in the observation space (with term Ri 2 ). Thus, the equation 29 
can be rewritten as : 

K = B ( Hf ) x 

Ri + HiBHf 

R. + B^BH^ J 1 (30) 

HiBHj + R i2 

H 2 BHf + 

To simplify the notation, let: 



Pi = HiBHf + Ri 
2 BH 2 r 



P 2 = H 2 BHo + R 



and: 



Pi 
P 2 



(31) 



(32) 



If, as assumed, the data assimilation procedure for each of the instruments 
has already been done, the terms Pi and P 2 are available as well are their 
inverse. Thus only extra diagonal term given by the following equation need to 
be treated: 

H 2 BHf + 6 

Also to simplify the notation, we denote: 

iT 



(33) 



S = H!BH^+R 12 (34) 

The other term of the matrix of Equation |33| can then be obtained by transposing 
the previous S matrix: 

S T = H 2 BHf + Rf 2 (35) 
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It worth to recall that B is a symmetric positive defined matrix, thus the equality 
B T = B can be used. 

Now, the aim is to put the part to invert in equation 30 under a specific 
form, that will ease the inversion. In this purpose, we write: 

U=(; t) (36) 

and: 

v=( s :;) 

If both matrix are multiplied in the following way, we notice a very useful 
property: 

uv =(s°- o) < 38 > 



In equation 38 we notice that extra diagonal terms can be expressed as a 
product of two matrix. Thus, rewriting the inverted term in equation |30[ we 
have to calculate the following: 

(P + UV)- 1 (39) 

This is from that expression that we realise inversion through Woodbury formula 
[25] written below: 

(P QR) 1 = P 1 + P X Q(I RP- 1 Q)- 1 RP- 1 (40) 

Changing sign on the left side, we obtain: 

(P + QR) 1 = P 1 P J Q(I + RP- 1 Q)- 1 RP- 1 (41) 

From equation |41[ we calculate the inversion needed in formula |39| First, we 
look after term P in equation |39| The P matrix is block-diagonal, thus we got : 

T £ ) 

Terms of the matrix are supposed to be known and stored a priori. Thus respect 
to formula 41 we got : 

(P + UV)- 1 = P 1 - P ^(I + VP _1 U) _1 VP- 1 (43) 

We take the following notation: 

C = P~ 1 U(I + VP^U^VP -1 (44) 

Thus a new formulation of that gain matrix K is obtained: 

K = BH(P 1 - C) (45) 

Within this equation, the term BHP 1 corresponds to join together the column 
of both matrix Ki and K2 from equation |25[ that are respectively of size n x d-y 
and n x d 2 . We will denote K this term BHP -1 , and thus we got: 

K = ( Kx K 2 ) = BHP 1 (46) 
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In the above equation, the index notice that the gain matrix is build on the 
basis of the matrix of each instruments set independently. Thus we can rewrite 



equation 45 in the following way : 

K = K BHC (47) 

The term BHC in equation [47] can be then interpreted as a correction to the 
matrix Ko, that is the simple coupling of the two original K matrix. Thus, we 
have build here a inversion method, where we can dissociate the contributions 
of each instrument and the one coming from their joint use. 



Looking more in detail to the matrix C, defined by equation 44 some points 
can be noticed. In the term (I — VP _1 U) _1 to be inverted, we recall that P 1 
matrix is of size (di + e? 2 ) x (di + d%). The V and U matrices are respectively 
of size Idi x [d\ + 1^2) and (di +c?2) x 2c?2. The matrix that need do be inverted 
in equation [44] is then of size 2di x 2di. If the size of the observation space 
associated to instrument 2 is smaller that the size of the observation space of 
instrument 1, it is far better to use this method. Of course the instrument 
denoted 1 and 2 can be ordered in such a way that di > d^, which ensure that 
the method always gives benefits. 

In the case of the permutation of instrument, we only know the analysis of 
the intermediate matrix for the all the instrument but one. Then the calculation 
of the analysis for only one instrument need to be performed before using this 
method. This correspond on overall to the inversion of an extra size matrix, 
which is rather cheap in calculation time. 
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